## Payson (2016) When Are Incumbents Held Accountable for Government Performance? 
## Evidence from U.S. School Districts
## Replication Code


## Load packages
library(ggplot2); library(plyr); library(effects); library(lmtest); library(RColorBrewer) 
library(plm); library(stargazer)

## Set working directory
setwd()

## Read in data
load("CaDistrictElections_2003-2012.RData")
str(dta)

## Format for panel analysis
pdta <- pdata.frame(dta, index = c("district", "year"))


############################################################################
## Tables (Main Paper)

## Table 1: Election descriptive statistics

elections <- ddply(dta, .(type), summarize,
							type = type[1], 
							num.race = round(length(unique(raceID)) / length(unique(year))),
							num.seat.mean = round(mean(num.seat, na.rm = T), 1),
							num.seat.sd = round(sd(num.seat, na.rm = T), 1),
							num.inc.mean = round(mean(num.inc, na.rm = T), 1),
							num.inc.sd = round(sd(num.inc, na.rm = T), 1),
							num.ch.mean = round(mean(num.ch, na.rm = T), 1),
							num.ch.sd = round(sd(num.ch, na.rm = T), 1),
							pct.inc.running.mean = round(mean(pct.inc.running, na.rm = T), 1),
							pct.inc.running.sd = round(sd(pct.inc.running, na.rm = T), 1),
							inc.voteshare.mean = round(mean(inc.voteshare, na.rm = T), 1),
							inc.voteshare.sd = round(sd(inc.voteshare, na.rm = T), 1),
							ch.voteshare.mean = round(mean(ch.voteshare, na.rm = T), 1),
							ch.voteshare.sd = round(sd(ch.voteshare, na.rm = T), 1),
							pct.inc.elected.mean = round(mean(pct.inc.elected, na.rm = T), 1),
							pct.inc.elected.sd = round(sd(pct.inc.elected, na.rm = T), 1),
							enroll.mean = round(mean(enroll, na.rm = T)),
							enroll.sd = round(sd(enroll, na.rm = T)),
							votes.mean = round(mean(num.voters, na.rm = T)),
							votes.sd = round(sd(num.voters, na.rm = T)),
							population.mean = round(mean(voting.pop, na.rm = T)),
							population.sd = round(sd(voting.pop, na.rm = T)),
							turnout.mean = round(mean(turnout, na.rm = T), 1),
							turnout.sd = round(sd(turnout, na.rm = T), 1),
							api.mean = round(mean(api.current, na.rm = T)),
							api.sd = round(sd(api.current, na.rm = T)),
							ayp.mean = round(mean(made.ayp1, na.rm = T), 2),
							ayp.sd = round(sd(made.ayp1, na.rm = T), 2))

## Table 1 (main paper)
stargazer(t(elections), summary = FALSE, type = "text",
					title = "California School Board Descriptive Statistics, 2003 - 2012")


####################################################################################	
## Table 2 Models: Changes in API score have a small and statistically insignificant 
## impact on the percentage of candidates who run

m1 <- plm(pct.inc.running ~ z.change + z.api + made.ayp + type + num.seat + year05 + 
							year07 + year08 + year09 + year10 + year11 + year12, 
							data = pdta, model = "within")
(coef1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC3", cluster = "group", adjust  = T)))

m2 <- plm(pct.inc.running ~ z.change*type + z.api + made.ayp + num.seat + year05 + 
							year07 + year08 + year09 + year10 + year11 + year12, 
							data = pdta, model = "within")
(coef2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef2)[14] <- "mid.interaction"
row.names(coef2)[15] <- "off.interaction"

m3 <- plm(pct.inc.running ~ api.change + api.current + made.ayp + type + num.seat + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
            data = pdta, model = "within")
(coef3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC3", cluster = "group", adjust  = T)))

m4 <- plm(pct.inc.running ~ api.change*type + api.current + made.ayp + num.seat + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
						data = pdta, model = "within")
(coef4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef4)[14] <- "mid.interaction"
row.names(coef4)[15] <- "off.interaction"


## Table 2 (main paper)
stargazer(coef3, coef4, coef1, coef2, type = "text",
						font.size = "footnotesize", label = "table:strategy",
						keep = c("api.change", "api.current", "z.change", "typeMidterm", "typeOff-Year",
											"mid.interaction", "off.interaction", "made.ayp", "z.api", "num.seat"),
						order = c("api.change", "z.change", "typeMidterm", "typeOff-Year",
											"mid.interaction", "off.interaction", "api.current",
											"z.api", "made.ayp", "num.seat"),
						title = paste("\\textbf{API Scores and Percent of Incumbents Seeking Re-Election.}",
												"Coefficients on both API changes and levels are small and imprecise,",
												"showing that district performance does not cause incumbents to",
												"strategically exit from school board races. The result holds if election",
												"types are considered separately."),												
						dep.var.labels = c("% Eligible Incumbents Seeking Re-Election (0 - 100)"),
						covariate.labels = c("One-Year API Change (Raw)", 
																"One-Year API Change (Standardized)",
																"Midterm Election", "Off-Year Election", 
																"API Change x Midterm", "API Change x Off-Year", 
																"API Level (Raw)", "API Level (Standardized)",
																"Made AYP", "Number of Open Seats"),
						add.lines = list(c("Year FE", "Y", "Y", "Y", "Y"), c("District FE", "Y", "Y", "Y", "Y"),
													c("", paste("Districts = ", length(fixef(m3))), 
															paste("Districts = ", length(fixef(m4))),
															paste("Districts = ", length(fixef(m1))),
															paste("Districts = ", length(fixef(m2)))),
													c("", paste("N = ", dim(m3$model)[1]), 
															paste("N = ", dim(m4$model)[1]),
															paste("N = ", dim(m1$model)[1]), 
															paste("N = ", dim(m2$model)[1]))))


##############################################################################
## Table 3 Models: Effect of API scores by election type -- interaction models

m1 <- plm(inc.voteshare ~ api.change*type + api.current + made.ayp + num.seat + 
						num.cand + year05 + year07 + year08 + year09 + year10 + year11 + year12, 
						data = pdta, model = "within")
(coef1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef1)[15] <- "mid.interaction"
row.names(coef1)[16] <- "off.interaction"

m2 <- plm(inc.voteshare ~ z.change*type + z.api + made.ayp + num.seat + num.cand + 
						year05 + year07 + year08 + year09 + year10 + year11 + year12, 
						data = pdta, model = "within")
(coef2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef2)[15] <- "mid.interaction"
row.names(coef2)[16] <- "off.interaction"

m3 <- plm(pct.inc.elected ~ api.change*type + api.current + made.ayp + num.seat + 
						num.cand + year05 + year07 + year08 + year09 + year10 + year11 + year12, 
						data = pdta, model = "within")
(coef3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef3)[15] <- "mid.interaction"
row.names(coef3)[16] <- "off.interaction"

m4 <- plm(pct.inc.elected ~ z.change*type + z.api + made.ayp + num.seat + num.cand + 
						year05 + year07 + year08 + year09 + year10 + year11 + year12, 
						data = pdta, model = "within")
(coef4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef4)[15] <- "mid.interaction"
row.names(coef4)[16] <- "off.interaction"


## Table 3 (Main Paper)
stargazer(coef1, coef2, coef3, coef4, type = "text",
						keep = c("api.change", "api.current", "z.change", "typeMidterm", "typeOff-Year",
											"mid.interaction", "off.interaction", "made.ayp", "z.api", "num.seat"),
						order = c("api.change", "z.change", "typeMidterm", "typeOff-Year",
											"mid.interaction", "off.interaction", "api.current",
											"z.api", "made.ayp", "num.seat"),
						title = paste("\\textbf{One-Year Change in Test Scores and Incumbent Performance.}",
												"In presidential election years, both raw and standardized one-year",
												"changes in API scores have strong, positive effects on incumbent",
												"voteshare and the percentage of incumbents that are re-elected.",
												"This relationship is almost completely flat in midterm years and even",
												"slightly negative in off-years"),											
						dep.var.labels = c("Avg. Incumbent Vote Share",
															"\\% Incumbents Re-Elected (0 - 100)"),
						covariate.labels = c("One-Year API Change (Raw)", 
															"One-Year API Change (Standardized)", "Midterm Election", 
															"Off-Year Election", "API Change x Midterm",
															"API Change x Off-Year", "API Level (Raw)", 
															"API Level (Standardized)", "Made AYP", 
															"Number of Open Seats"),
						add.lines = list(c("Year FE", "Y", "Y", "Y", "Y"), c("District FE", "Y", "Y", "Y", "Y"),
													c("", paste("Districts = ", length(fixef(m1))), 
															paste("Districts = ", length(fixef(m2))),
															paste("Districts = ", length(fixef(m3))),
															paste("Districts = ", length(fixef(m4)))),
													c("", paste("N = ", dim(m1$model)[1]), 
															paste("N = ", dim(m2$model)[1]),
															paste("N = ", dim(m3$model)[1]), 
															paste("N = ", dim(m4$model)[1]))))


###########################################################################
## Table 4 Models: Voter turnout is lower in midterm and off-year elections

m <- plm(turnout ~ type, data = pdta, model = "pooling")
summary(m)
(coef1 <- coeftest(m, vcov = vcovHC(m, type = "HC3", cluster = "group", adjust  = T)))

m1 <- plm(turnout ~ type, data = pdta, model = "within")
summary(m1)
(coef2 <- coeftest(m1, vcov = vcovHC(m1, type = "HC3", cluster = "group", adjust  = T)))


## Table 4 (main paper)
stargazer(coef1, coef2, type = "text",
						title = "Turnout By Election Type",
						dep.var.labels = c("Voting Age Turnout"),
						column.labels = c("Pooled", "Fixed Effects"),
						covariate.labels = c("Midterm", "Off-Year", "(Intercept)"),
						add.lines = list(c("District FE", "", "Y"),
													c("", "", paste("Districts = ", length(fixef(m1)))),
													c("", paste("N = ", dim(m$model)[1]),
															paste("N = ", dim(m1$model)[1]))))

						
############################################################################
## Appendix Tables													
############################################################################									

## Table A1 Models: Poisson regression

m1 <- glm(num.inc ~ z.change + z.api + num.seat + year + district, 
						data = dta, family = "poisson")
summary(m1)$coefficients[1:4, ]

m2 <- glm(num.inc ~ z.change*type + z.api + num.seat + year + district, 
						data = dta, family = "poisson")
summary(m2)$coefficients[c(1:6, 849:850), ] 


## Table A1 (Appendix)
stargazer(m1, m2, type = "text",
					keep = c("z.change", "type", "z.change*type", "z.api", "num.seat"))


##################################################################################
## Table A2 Appendix Models: no effect of API on incumbent voteshare when election
## types are pooled 

m1 <- plm(inc.voteshare ~ api.change + api.current + made.ayp + num.cand + num.seat + year, 
							data = pdta, model = "within")
(coef1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC3", cluster = "group", adjust  = T)))

m2 <- plm(inc.voteshare ~ z.change + z.api + made.ayp +num.cand + num.seat + year, 
							data = pdta, model = "within")
(coef2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC3", cluster = "group", adjust  = T)))

m3 <- plm(pct.inc.elected ~ api.change + api.current + made.ayp + num.cand + num.seat + year, 
								data = pdta, model = "within")
(coef3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC3", cluster = "group", adjust  = T)))

m4 <- plm(pct.inc.elected ~ z.change + z.api + made.ayp + num.cand + num.seat + year, 
							data = pdta, model = "within")
(coef4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC3", cluster = "group", adjust  = T)))


## Table A2 (Appendix)
stargazer(coef1, coef2, coef3, coef4, type = "text",
						keep = c("z.change", "z.api", "api.change", "api.current", 
											"made.ayp", "num.cand", "num.seat"),
						order = c("api.change", "api.current", "z.change", "z.api", 
											"made.ayp", "num.cand", "num.seat"),				
						title = "One-Year Change in Test Scores and Incumbent Performance (All Years)",
						column.labels = c("Avg. Incumbent Vote Share (0 - 100)", 	
															"\\% Incumbents Re-Elected (0 - 100)"),
						column.separate = c(2, 2),
						covariate.labels = c("One-Year Change in API Change (Raw)", "API Level (Raw)",
															"One-Year API Change (Standardized)", "API Level (Standardized)",
															"Made AYP", "Number of Candidates", "Number of Seats"),
						add.lines = list(c("District FE", "Y", "Y", "Y", "Y"), c("Year FE", "Y", "Y", "Y", "Y"),
													c("", paste("Districts = ", length(fixef(m1))), 
															paste("Districts = ", length(fixef(m2))), 
															paste("Districts = ", length(fixef(m3))), 
															paste("Districts = ", length(fixef(m4)))),
													c("", paste("N = ", dim(m1$model)[1]), 
															paste("N = ", dim(m2$model)[1]), 
															paste("N = ", dim(m3$model)[1]), 
															paste("N = ", dim(m4$model)[1]))))


############################################################################
## Table A3 Models: One-Year Change in Test Scores and Incumbent Performance
## (Additional Covariates)

m1 <- plm(pct.inc.elected ~ z.change*type + z.api + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
            data = pdta, model = "within")
(coef1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef1)[12] <- "mid.interaction"
row.names(coef1)[13] <- "off.interaction"

m2 <- plm(pct.inc.elected ~ z.change*type + z.api + made.ayp + num.seat + num.cand + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
            data = pdta, model = "within")
(coef2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef2)[15] <- "mid.interaction"
row.names(coef2)[16] <- "off.interaction"

m3 <- plm(pct.inc.elected ~ z.change*type + z.api + made.ayp + 
            log.enroll + pct.low.ses + pct.minority + num.seat + num.cand + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
            data = pdta, model = "within")
(coef3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef3)[18] <- "mid.interaction"
row.names(coef3)[19] <- "off.interaction"


## Table A3 (Appendix)
stargazer(coef1, coef2, coef3, type = "text", 
						keep = c("z.change", "typeMidterm", "typeOff-Year", "log.enroll", "pct.low.ses",
											"mid.interaction", "off.interaction", "made.ayp", "z.api", "num.seat",
											"num.cand", "pct.minority"),
						order = c("z.change", "typeMidterm", "typeOff-Year",
											"mid.interaction", "off.interaction", "z.api", "made.ayp", "num.seat",
											"num.cand", "log.enroll", "pct.low.ses", "pct.minority"),
						title = paste("\\textbf{One-Year Change in Test Scores and Incumbent Performance.}",
												"In presidential election years, both raw and standardized one-year",
												"changes in API scores have strong, positive effects on incumbent",
												"voteshare and the percentage of incumbents that are re-elected.",
												"This relationship is almost completely flat in midterm years and even",
												"slightly negative in off-years"),											
						dep.var.labels = c("\\% Incumbents Re-Elected (0 - 100)"),
						covariate.labels = c("One-Year API Change (Standardized)",
																"Midterm Election", "Off-Year Election", "API Change x Midterm",
																"API Change x Off-Year", "API Level (Standardized)",							
																"Made AYP", "Number of Seats", "Numer of Candidates",
																"District Enrollment (Log)", "District \\% Low SES", 
																"District \\% Minority"),
						add.lines = list(c("Year FE", "Y", "Y", "Y"), c("District FE", "Y", "Y", "Y"),
													c("", paste("Districts = ", length(fixef(m1))), 
															paste("Districts = ", length(fixef(m2))),
															paste("Districts = ", length(fixef(m3)))),
													c("", paste("N = ", dim(m1$model)[1]), 
															paste("N = ", dim(m2$model)[1]),
															paste("N = ", dim(m3$model)[1]))))


###########################################
## Table A4 Models: Random District Effects

r1 <- plm(inc.voteshare ~ z.change*type + z.api + num.seat + num.cand + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
          data = pdta, model = "random")
(rcoef1 <- coeftest(r1, vcov = vcovHC(r1, type = "HC0", cluster = "group", adjust  = T)))

r2 <- plm(pct.inc.elected ~ z.change*type + z.api + num.seat + num.cand + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
          data = pdta, model = "random")
(rcoef2 <- coeftest(r2, vcov = vcovHC(r2, type = "HC0", cluster = "group", adjust  = T)))


## Table A4 (Appendix)
stargazer(rcoef1, rcoef2, type = "text",
						title = "One-Year Change in Test Scores and Incumbent Performance: 
										Random District Effects",
						keep = c("z.change", "z.change:typeMidterm", "z.change:typeOff-Year",
											"typeMidterm", "typeOff-Year", "z.api", "num.seat", "num.cand"),
						order = c("z.change", "z.change:typeMidterm", "z.change:typeOff-Year",
											"typeMidterm", "typeOff-Year", "z.api", "num.cand", "num.seat"),
						column.labels = c("Avg. Incumbent Vote Share", "% Incumbents Re-Elected", 
															"All Incumbents Re-Elected"),
						covariate.labels = c("One-Year API Change", "API Change x Midterm Election", 
															"API Change x Off-Year Elections", "Midterm Election",
															"Off-Year Election", "API Level (Standardized)",
															"Number of Candidates", "Number of Seats"))


########################################################
## Table A5 Models: API Levels and Incumbent Performance

m1 <- plm(inc.voteshare ~ api.current*type + api.change + made.ayp + num.seat + 
						num.cand + year05 + year07 + year08 + year09 + year10 + year11 + year12, 
						data = pdta, model = "within")
(coef1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef1)[15] <- "mid.interaction"
row.names(coef1)[16] <- "off.interaction"

m2 <- plm(inc.voteshare ~ z.api*type + z.change + made.ayp + num.seat + num.cand + 
						year05 + year07 + year08 + year09 + year10 + year11 + year12, 
						data = pdta, model = "within")
(coef2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef2)[15] <- "mid.interaction"
row.names(coef2)[16] <- "off.interaction"

m3 <- plm(pct.inc.elected ~ api.current*type + api.change + made.ayp + num.seat + 
						num.cand + year05 + year07 + year08 + year09 + year10 + year11 + year12, 
						data = pdta, model = "within")
(coef3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef3)[15] <- "mid.interaction"
row.names(coef3)[16] <- "off.interaction"

m4 <- plm(pct.inc.elected ~ z.api*type + z.change + made.ayp + num.seat + num.cand + 
						year05 + year07 + year08 + year09 + year10 + year11 + year12, 
						data = pdta, model = "within")
(coef4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC3", cluster = "group", adjust  = T)))
row.names(coef4)[15] <- "mid.interaction"
row.names(coef4)[16] <- "off.interaction"


## Table A5 (Appendix)
stargazer(coef1, coef2, coef3, coef4, type = "text",
						keep = c("api.change", "api.current", "z.change", "typeMidterm", "typeOff-Year",
											"mid.interaction", "off.interaction", "made.ayp", "z.api", "num.seat"),
						order = c("api.current", "z.api", "typeMidterm", "typeOff-Year",
											"mid.interaction", "off.interaction", "z.change",
											"api.change", "made.ayp", "num.seat"),
						title = paste("While one-year changes in API have a strong and consistent effect on",	
												"incumbent performance depending on when an election is held,",
												"API levels have no effect across election types."),										
						dep.var.labels = c("Avg. Incumbent Vote Share",
															"\\% Incumbents Re-Elected (0 - 100)"),
						covariate.labels = c("API Level (Raw)", "API Level (Standardized)",
																"Midterm Election", "Off-Year Election", "API Level x Midterm",
																"API Level x Off-Year", "One-Year API Change (Raw)", 
																"One-Year API Change (Standardized)",
																"Made AYP", "Number of Open Seats"),
						add.lines = list(c("Year FE", "Y", "Y", "Y", "Y"), 
						                 c("District FE", "Y", "Y", "Y", "Y"),
													c("", paste("Districts = ", length(fixef(m1))), 
															paste("Districts = ", length(fixef(m2))),
															paste("Districts = ", length(fixef(m3))),
															paste("Districts = ", length(fixef(m4)))),
													c("", paste("N = ", dim(m1$model)[1]), 
															paste("N = ", dim(m2$model)[1]),
															paste("N = ", dim(m3$model)[1]), 
															paste("N = ", dim(m4$model)[1]))))
															
															
#########################################################################
## Figures (Main Paper)															
#########################################################################

## Read in data for Figures 1 & 2
dist <- read.csv("CaDistrictPerformance_2003-2012.csv")

## Figure 1 setup
mean.api <- tapply(dist$api.prior, dist$year, function (x) mean(x, na.rm = T))
sd.api <- tapply(dist$api.prior, dist$year, function (x) sd(x, na.rm = T))
n <- unname(table(dist$year))

lower <- as.numeric(mean.api - (1.96)*(sd.api / sqrt(n)))
upper <- as.numeric(mean.api + (1.96)*(sd.api / sqrt(n)))

figdta <- data.frame(year = seq(2003, 2012), mean.api, lower, upper)

## Figure 1 (Main Paper)
p <- ggplot(figdta, aes(factor(year), mean.api, group = 1))
p + geom_point() +
		geom_line() +
		geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.05) +
		theme_bw(base_size = 12, base_family = "serif") +
		labs(y = "Average API Score", x = "Year")
		

## Figure 2 setup
figdta$ackerman <- dist$api.current[dist$district == "Ackerman Elementary"]
figdta$alpaugh <- dist$api.current[dist$district == "Alpaugh Unified"]

## Figure 2 (main paper)
p <- ggplot(figdta, aes(factor(year), mean.api, group = 1))
p + geom_point() +
		geom_point(aes(x = 2, y = figdta$ackerman[2]), shape = 1, size = 5) +
		geom_point(aes(x = 10, y = figdta$ackerman[10]), shape = 1, size = 5) +
		geom_point(aes(x = 4, y = figdta$alpaugh[4]), shape = 1, size = 5) +
		geom_point(aes(x = 10, y = figdta$alpaugh[10]), shape = 1, size = 5) +
		geom_line() +
		geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.05) +
		geom_line(aes(x = factor(year), y = ackerman)) +
		geom_line(aes(x = factor(year), y = alpaugh)) +
		geom_segment(aes(x = 3.5, xend = 3.5, y = 510, yend = 610), linetype = 2) +
		geom_segment(aes(x = 9.5, xend = 9.5, y = 650, yend = 750), linetype = 2) +
		geom_segment(aes(x = 1.5, xend = 1.5, y = 775, yend = 875), linetype = 2) +
		geom_segment(aes(x = 9.5, xend = 9.5, y = 825, yend = 925), linetype = 2) +
		theme_bw(base_size = 12, base_family = "serif") +
		labs(y = "API Score", x = "Year") +
		annotate("text", x = 6, y = 850, label = "Ackerman Elementary", size = 4, family = "serif") +
		annotate("text", x = 7, y = 630, label = "Alpaugh Unified", size = 4, family = "serif") +
		annotate("text", x = 2.8, y = 740, label = "All District Average", size = 4, family = "serif") +
		geom_point(aes(x=7.60, y=550), shape = 1, size = 5) +
		annotate("text", x = 8.67, y = 550, label = "= Election Year", size = 4, family = "serif") 


##################################################################
## Figure 3 (Top): API change and incumbent vote share, all years 
ggplot(pdta, aes(z.change, inc.voteshare)) +
	geom_jitter(alpha = .3) +
	ylim(0, 100) +
	scale_x_continuous(limits = c(-5, 5), breaks = c(-4, -2, 0, 2, 4)) +
	geom_smooth(method = "lm", colour = "black") +	
	theme_bw(base_family = "serif") +
	labs(y = "Average Incumbent Vote Share", x = "One Year Change in API Score (Standardized)") +
	theme(strip.background = element_rect(fill = "white"))

## Figure 3 (Bottom): API change and incumbent vote share, subset by election type
ggplot(pdta, aes(z.change, inc.voteshare)) +
	geom_point(alpha = .3) +
	facet_wrap(~type) +
	ylim(0, 100) +
	scale_x_continuous(limits = c(-5, 5), breaks = c(-4, -2, 0, 2, 4)) +
	geom_smooth(method = "lm", colour = "black") +	
	theme_bw(base_family = "serif") +
	labs(y = "Average Incumbent Vote Share", x = "One Year Change in API Score (Standardized)") +
	theme(strip.background = element_rect(fill = "white"))
	
	
##############################################
## Figure 4 models and setup: Marginal Effects

m <- lm(pct.inc.elected ~ z.change*type + z.api + num.seat + num.cand + 
          year05 + year07 + year08 + year09 + year10 + year11 + year12 + district, data = pdta)

plot <- data.frame(Type = c(rep("Presidential", 61), rep("Midterm", 61), rep("Off-Year", 61)),
									Lower = c(summary(effect("z.change*type", m, 
															xlevels = list(z.change = seq(-3, 3, .1))))$lower[,1],
														summary(effect("z.change*type", m, 
															xlevels = list(z.change = seq(-3, 3, .1))))$lower[,2],
														summary(effect("z.change*type", m, 
															xlevels = list(z.change = seq(-3, 3, .1))))$lower[,3]),
								 	Upper = c(summary(effect("z.change*type", m, 
								 							xlevels = list(z.change = seq(-3, 3, .1))))$upper[,1],
														summary(effect("z.change*type", m, 
															xlevels = list(z.change = seq(-3, 3, .1))))$upper[,2],
														summary(effect("z.change*type", m, 
															xlevels = list(z.change = seq(-3, 3, .1))))$upper[,3]),
									Effect = c(summary(effect("z.change*type", m, 
										xlevels = list(z.change = seq(-3, 3, .1))))$effect[,1],
														summary(effect("z.change*type", m, 
															xlevels = list(z.change = seq(-3, 3, .1))))$effect[,2],
														summary(effect("z.change*type", m, 
															xlevels = list(z.change = seq(-3, 3, .1))))$effect[,3]),
									API = rep(seq(-3, 3, .1), 3))

plot$Type <- factor(plot$Type, levels = c("Presidential", "Midterm", "Off-Year"))


## Figure 4: Marginal Effect of Test Score Changes on Incumbent Re-election
ggplot(plot, aes(API, Effect)) +
	geom_line(aes(x = API, y = Effect)) +
	geom_line(aes(x = API, y = Lower), linetype = 2) +
	geom_line(aes(x = API, y = Upper), linetype = 2) +
	facet_wrap(~Type) +
	ylim(20, 100) +
	theme_bw(base_family = "serif") +
	theme(plot.background = element_blank(),
   				panel.grid.major = element_blank(),
   				panel.grid.minor = element_blank(),
   				strip.background = element_rect(fill = "white")) +
	labs(y = "% Incumbents Re-Elected", x = "API Score (Standardized)")				


#####################################################
## Figure 5 models: Fixed vs. Random Coefficient Plot

m1 <- plm(inc.voteshare ~ z.change*type + z.api + num.seat + num.cand + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
          data = pdta, model = "within")
(coef1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC0", cluster = "group", adjust  = T)))

m2 <- plm(pct.inc.elected ~ z.change*type + z.api + num.seat + num.cand + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
          data = pdta, model = "within")
(coef2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC0", cluster = "group", adjust  = T)))

r1 <- plm(inc.voteshare ~ z.change*type + z.api + num.seat + num.cand + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
          data = pdta, model = "random")
(rcoef1 <- coeftest(r1, vcov = vcovHC(r1, type = "HC0", cluster = "group", adjust  = T)))

r2 <- plm(pct.inc.elected ~ z.change*type + z.api + num.seat + num.cand + 
            year05 + year07 + year08 + year09 + year10 + year11 + year12, 
          data = pdta, model = "random")
(rcoef2 <- coeftest(r2, vcov = vcovHC(r2, type = "HC0", cluster = "group", adjust  = T)))

type <- c(rep(c("Fixed", "Random"), 6))

beta <- c(coef1[1,1], rcoef1[2,1], coef1[14,1], rcoef1[15, 1], coef1[15,1], rcoef1[16,1],
						coef2[1,1], rcoef2[2,1], coef2[14,1], rcoef2[15, 1], coef2[15,1], rcoef2[16,1])

se <- c(coef1[1,2], rcoef1[2,2], coef1[14,2], rcoef1[15, 2], coef1[15,2], rcoef1[16,2],
					coef2[1,2], rcoef2[2,2], coef2[14,2], rcoef2[15, 2], coef2[15,2], rcoef2[16,2])

y <- c(23, 22, 19, 18, 15, 14, 10, 9, 6, 5, 2, 1)

figdata <- data.frame(type, beta, se, y)


## Figure 5: Fixed vs. Random District Effects
ggplot(figdata, aes(x = beta, y = y, lty = type)) + 
	geom_point(size = 2) +
	scale_linetype_discrete(name = "Model Type", 
												breaks = c("Fixed", "Random"), 
												labels = c("Fixed Effects", "Random Effects")) +
	geom_errorbarh(aes(xmin = beta - (1.96 * se), xmax = beta + (1.96 * se)), 
											size = .5, height = .8) +
	geom_vline(xintercept = 0, linetype = 2) +
	geom_hline(yintercept = 12, linetype = 1, size = .25) +
	scale_y_continuous(breaks = c(1.5, 5.5, 9.5, 14.5, 18.5, 22.5), 
											labels = rep(c("One Year Change in API x Off-Year Elections", 
																		"One Year Change in API x Midterm Elections", 
																		"One Year Change in API x Presidential Elections"), 2)) + 
	theme_bw(base_family = "serif") +  labs(y = " ", x = "") +
	theme(plot.background = element_blank(),
   				panel.grid.major = element_blank(),
   				panel.grid.minor = element_blank()) +
	annotate("text", x = -7, y = 24, label = "Incumbent Vote Share", size = 3, family = "serif") +
	annotate("text", x = -5.5, y = 11, label = "Percent of Incumbents Re-Elected", 
						size = 3, family = "serif")


#######################
## Figure A1 (Appendix)

ggplot(dta, aes(z.change, num.inc)) +
	geom_jitter(alpha = .3) + 
	stat_smooth(method = "lm", colour = "black") +
	xlim(-2, 2) +
	theme_bw(base_family = "serif") +
	labs(y = "Number of Incumbents Running", x = "One Year Change in API Score (Standardized)")